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Abstract 


Moran's  I  statistic  has  been  used  by  ecologists  and  geographers 
alike  to  test  for  the  presence  of  spatial  autocorrelation  in  a  single 
variable  over  a  two-dimensional  plane.  In  this  paper,  we  provide  the 
derivation  of  the  expected  value  and  variance  of  a  bivariate  version 
of  Moran's  /  for  use  with  multivariate  data  under  the  assumption  of 
spatial  independence.  We  also  demonstrate  that  Moran's  univariate 
I  statistic  is  a  special  case  of  Moran's  bivariate  IYZ.  Results  of  an 
extensive  Monte  Carlo  study  show  that  the  expected  value  and 
variance  are  reliable  for  several  data  sets  with  moderate  sample  sizes 
[n  =  40  and  127)  and  varying  degrees  of  correlation  among  different 
bivariate  surfaces.  For  small  sample  sizes  [n  <  40)  at  least  5,000 
simulations  should  be  used  to  estimate  the  P-value  to  test  the  null 
hypothesis  of  no  spatial  autocorrelation  between  two  variables. 
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Expected  Value  and  Variance  of  Moran's  Bivariate 
Spatial  Autocorrelation  Statistic  for  a  Permutation  Test 


Raymond  L.  Czaplewski  and  Robin  M.  Reich 


Introduction 

In  analyzing  data  over  large  geographical  areas,  one 
is  often  interested  in  knowing  whether  or  not  there  is  a 
spatial  relationship  among  the  data,  and  if  there  is,  can 
one  obtain  information  on  the  process  responsible  for 
generating  the  observed  pattern.  Geostatistical  tech- 
niques, such  as  Moran's  statistic  of  spatial  autocorrelation 
(Cliff  and  Ord  1973),  spectral  analysis  (Ripley  1981), 
and  the  Mantel  test  (Mantel  1967),  are  just  a  few  proce- 
dures which  permit  tests  of  hypotheses  regarding  spa- 
tial patterns.  These  tests  not  only  describe  the  structure 
of  the  spatial  pattern  but  they  are  also  capable  of  detect- 
ing the  presence  of  directional  components  at  various 
scales  (Legendre  and  Fortin  1989).  Although  these 
techniques  have  been  applied  to  univariate  data,  they 
can  be  extended  to  handle  multivariate  data  as  well. 

In  this  paper  we  provide  the  derivation  of  a  bivariate 
version  Moran's  /,  denoted  IYZ,  for  evaluating  spatial 
autocorrelation  under  the  null  hypothesis  that  the  mag- 
nitude of  spatial  autocorrelation  is  no  greater  than  that 
expected  by  the  chance  distribution  over  space  of  inde- 
pendent events.  We  also  demonstrate  that  Moran's 
univariate  /statistic  (Cliff  and  Ord  1973)  is  a  special  case 
of  Moran's  bivariate  7YZ.  Finally,  we  discuss  the  statis- 
tical properties  of  Moran's  bivariate  statistic  using  data 
from  the  USDA  Forest  Service,  Southeastern  Forest 
Experiment  Station  Inventory  and  Analysis  Program, 
for  natural,  undisturbed  shortleaf  pine  stands  [Pinus 
echinata)  in  northern  Georgia. 


Bivariate  Moran's  IYZ 

Define  statistic  IYZ  as  an  index  of  spatial 
autocorrelation  between  two  variables  ( Y and  Z) ,  each  of 
which  are  observed  at  n  locations  (IYZ  is  analogous  to 
Moran's  I  for  a  univariate  variable): 

n  n 

iYZ  =  1=1  ]=i 

'**  ,  [1] 

where 

w..  =     a  scalar  that  quantifies  the  degree  of 
spatial  association  or  proximity  between 


locations  i  and  /  (e.g.,  inverse  distance 
between  locations  i  and  /,  or  a  0-1  vari- 
able indicating  that  locations  i  and  /  are 
within  some  distance  range  of  each 
other);  w..  does  not  necessarily  equal  w..; 
w..  -  0  for  all  i 
Y.  =  the  observed  value  of  variable  Y  for  plot 
i,  transformed  so  its  mean  is  zero 

n 


Z.  =      the  observed  value  of  variable  Z  for  plot 
/,  transformed  so  its  mean  is  zero 

n 

I  z        n     _  n 

IT-     p=1  P 

n  n 

W=       II  w.. 
1=1  ,=i 

n 

I  y2 

mv2=     p=1  P  =Var(Y) 
n 

n 

I  z2 

m72=     p=1  P     =  Var(Z) 
L  n 

i=  (1,  2, n) 
/  =        (1,  2,  ...,n) 

n  =       number  of  observed  locations,  n  >  4. 

The  denominator  in  Eq.  1  makes  Jyz  a  dimensionless 
statistic  that  can  be  interpreted  as  a  weighted  correla- 
tion coefficient  in  space  between  variables  Y  and  Z. 
Thus,  one  would  expect  IYZ  to  range  over  the  interval  of 
-1  to  1.  Large  positive  values  of  Iyz  would  indicate  a 
positive  spatial  autocorrelation  between  the  two  re- 
sponse surfaces,  while  a  large  negative  value  would 
indicate  a  negative  spatial  autocorrelation.  Cliff  and 
Ord  (1981)  point  out  that  in  general,  the  upper  bound  for 
Moran's  univariate  lis  usually  less  than  1,  although  it 
can  exceed  1  for  an  irregular  pattern  of  weights,  w..,  if 
extreme  values  are  heavily  weighted.  Thus,  depending 
on  the  weights  used,  it  is  possible  for  Moran's  bivariate 
/V7  to  exceed  the  limits  of  -1  and  1. 
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Assume  the  bivariate  observation  {y.,  z.}  for  location 
i  is  a  random,  spatially  independent  drawing  from  one 
(or  separate  identical)  population(s),  and  the  joint  dis- 
tribution function(s)  for  Y and  Z  are  unknown.  The  null 
hypothesis  is  that  there  is  no  spatial  autocorrelation, 
without  any  parametric  assumptions  regarding  the  joint 
distribution  of  Y and  Z.  Under  this  null  hypothesis,  we 
may  consider  the  set  of  all  nl  random  permutations,  in 
which  the  bivariate  observations  (yk,  zk},k=  {1,  2,  n}, 
are  randomly  assigned  without  replacement  to  location 
i,  i  =  {1,  2, n],  regardless  of  the  original  spatial  loca- 
tion of  the  kth.  observation.  Under  the  null  hypothesis 
of  no  spatial  autocorrelation,  each  bivariate  observation 
is  equally  likely  at  any  observation  location,  i.e.,  the 
random  assignment  of  any  bivariate  observation  {yk,  zk} 
to  location  i  is  equally  likely: 

prob(location  i  has  observation  {yk,  xk})  =  1/n, 
for  all  i  and  k. 

In  the  following  sections,  we  derive  the  exact  mean 
and  variance  of  Moran's  bivariate  JYZ  under  the  null 
hypothesis  for  the  n\  random  permutations,  and  condi- 
tional upon  the  bivariate  observations  at  the  n  locations. 
The  derivations  closely  parallel  those  for  Moran's 
univariate  /  given  by  Cliff  and  Ord  (1973)  under  their 
randomization  assumption.  The  derivation  proceeds  in 
the  next  sections  as  follows:  various  moments  of  the 
joint  distributions  of  Y  and  Z  are  derived;  these  expec- 
tations are  used  to  compute  the  expected  values  of  JYZ 
and  1^  (Eqs.  25  and  30);  simplifications  are  introduced 
to  make  computations  of  the  expected  value  of  more 
efficient  and  less  prone  to  numerical  problems  (Eq.  58), 
and  then  produce  an  exact  equation  for  the  variance  of 
JYZ  (Eq.  60);  and  we  then  show  that  the  mean  and 
variance  of  Moran's  univariate  7  is  a  special  case  of 
Moran's  bivariate  Iyz  which  we  present  here  (Eqs.  66,  68, 
and  70).  The  expected  value  is  denoted  by  E[«];  the 
mean  is  the  expected  value  of  Jyz,  i.e.,  E[/yz];  and  the 
variance  is  denoted  Var(/YZ)  and  is  equal  to  E [7^1  -  E[Jyz]2 
by  definition. 

Moments  of  Joint  Probability  Density  Function 

Let  Y.  be  the  value  of  the  y-  variable  for  location  i,  and 
Z.  be  the  value  of  the  z-variable  for  location  i  during  any 
one  of  the  nl  random  permutations.  Under  the  null 
hypothesis,  the  expected  value  of  the  variable  Y.,  taken 
over  all  possible  n\  permutations,  is  by  definition: 

E[Y)  =  7l  P(r=7l)  +  y2  P(l>y2)  +  ...  +  yn  P(l>yn)  [2] 

where  P[Y.=y  )  denotes  the  probability  that  the 
permutated  random  variable  Y.,  i<i<n,  equals  the 
realized  observation  y ,  i<  p  <  n,  over  all  possible  nl 
permutations.  There  are  (n  -  1)!  permutations  among  all 
possible  n!  permutations  in  which  Y.  =  y  ,[n  -  1)\  per- 


mutations in  which  Y.  =  y2,  and  so  forth.  Therefore, 
P(y.=yp)  =  [n  -  l)\ln\  =  1/n,  and  Eq.  2  simplifies  to: 


E[Y.]  =  I  yP 

1  p=l 


(n-1) 


n  ! 


n 

I  yF 

n/  p=1 


[3] 


Since  the  y  values  are  transformed  to  center  its  mean  on 
zero  over  the  n  observations  without  loss  of  generaliza- 
tion, we  can  define  the  E[Y.]  in  Eq.  1,  to  be  equal  to: 


n 

E  y 

E[Y.]  =  =  0 

n 


[4] 


Likewise,  the  expected  values  of  the  variables  Y.  and  Z. 
and  their  cross-products,  taken  over  all  possible  nl 
permutations,  are: 


E[Z.]  =  P=i  p  =0 


n 


n 

E[Y.Z.]  =      YpZp  =  Cov(rZ)  =  m 
n 

n 

E[Y?]  =  hy*  =Var(Y)  =  mY5 
n 


y  z2 

E[Zf]=  P=i  P    =Var(Z)  =  m2 
n 


E[YfZf]=hl^L  =my2z2 
n 


YZ 


[5] 


[6] 


[7] 


[8] 


[9] 


for  all  i. 

In  any  random  permutation  under  the  null  hypoth- 
esis, a  particular  location  i  can  have  value  Y.=yk  with 
probability  and  any  neighboring  unit /  (;W)  can  take 
any  value  Z.=zm  other  than  Z.=zk  with  probability 
l/(n-l);  therefore: 

n 

E[y.z.]  =  iZL  fz1  +  ---  +  zi,1  +  zi+1  +  ...zn\ 


n 


X  I  Yi  Zj 
=  i=i  i=l 


(n-1) 


n(n-l) 


i=l 

[10] 


n(n-l) 


The  sum  of  all  zp's  in  Eq.  10  equals  zero  by  definition 
(Eqs.  1  and  5),  and 


2 


e[y.z.]  =-x 

1  I  i=l 


n(n-l) 


[11] 


n  n 

Since  X  Yft  in  Eq.  11  equals  1^  ypzp  in  Eq.  6, 


1=1 


E[y.z]  =  "nmvz  -  'mYZ 
n(n-l)  n-1 


[12] 


where  mYZ  is  defined  in  Eq.  6.  Eq.  12  applies  to  all 
possible  pairs  of  locations  i  and  /,  i-tj,  in  the  nl  random 
permutations. 

In  a  similar  fashion,  EfY^Z2]  can  be  derived  for  fej 
using  Eqs.  7,  8,  and  9: 


n  n 

X  X  y2z2 


n(n-l) 


n 

i=l 


X  zp  N 


p=l 


n(n-l) 

n  n 

=  X  y?  (nmz2^  -X  y?  H 


i=l 


i=l 


n(n-l) 

_  nmy2  (nmzi)-nmY2z2  _  nmy2  mz2  -my2z2  [\2>\ 
n(n-l)  n-1 

where  mY2,  mz2,  and  mY2Z2  are  defined  in  Eqs.  7,  8,  and  9 
respectively.  Next,  the  E[Y.ZY.Z],  for  Utj,  can  be  de- 
rived in  a  fashion  similar  to  Eq.  13  using  Eqs.  6  and  9: 


n  n 

E[Y.Z.YZ]  =  X  X  y^ty^) 
1   )    i   1       i=l  j=l 

N  

n(n-l) 


n 

X  y>z. 

i=l 


n(n-l) 

£  y.z.(nmYZ)-X  yf  [z|] 

i  =  l 


i=l 


n(n-l) 

nmYZ  (nmYZ )  -  nmyz  z2_  nmyz  -  mY2z2 


n(n-l) 


[14] 


n-1 


where  mYZ  and  mY2z2are  defined  in  Eqs.  6  and  9,  respec- 
tively. 

Similarly,  the  E[  Y!ZZJ,for  itj,  &k,  and  jtk,  is  derived  as 
follows  using  Eqs.  7,8,  and  9: 


n    n  n 

XXX  It*?* 
E[Y!Z.Zk]=i=i  j=ik=i 

1    >    K  i*j*k 


n(n-l)  (n-2) 


n  n 

XXyfz 

i=i  j=i 


n(n-l)  (n-2) 


n    n  n  n 

-X  X  yftzj  -  X  X  y2*2 

i=l  j=l  i=l  j  =  l 

i*j  if)  

n(n-l)  (n-2) 


n 

-Xyfzi 

i=l 

=  1=1 


n(n-l)  (n-2) 

n 

X  yfz.[-z.]  -  X  yf[nmz2  -  z2] 


i=l 


n(n-l)  (n-2) 

n  n  n 

X  v2z2-nm72  X  v2+X  y2z2 
=  i=iJi  1      z  i=i  1   i=r  1  1 


n(n-l)  (n-2) 

nmY2z2  -  nmZ2  nmY2+  nmy2z2 
n(n-l)  (n-2) 

2my2z2  -  nmy2  mz2 


(n-1)  (n-2) 


[15] 


The  E[y.ymZ2],  for  i±m,  and  j*m,  is  derived  similar  to 
Eq.  15:  ' 


n    n  n 


xxx  yi  ymzi 


E[Y.Y  Z2]  =i=ij=i  m=i 

1    m   /  i*j*m 


n(n-l)  (n-2) 


2mY2Z2  -  nmy2  mz2 
(n-1)  (n-2) 


[16] 


E[Y.Z.Y.Zk],  for  tej,  i*ic,  and  fek,  is  derived  similar 
to  Eq.  15  using  Eqs.  6  and  9: 

n    n  n 

E[Y.Z.Y.ZV]  =  2  X   X  y^ 

1   i    i   k      i=l  j=l  k=l 

i*j*k  


n(n-l)  (n-2) 


n  n 
1=1 1=1 


n(n-l)  (n-2) 


n  n  n 

=  111 

i=l  j=l  m=l 
i*j*m 


I   Zn  -Z.-Z.-Z 


n(n-l)  (n-2)  (n-3) 


n  n 


n  n 


=  1=1  )=i       1=1  j=i 


n(n-l)  (n-2) 


n 

1=1 


n 


y.zz 


n(n-l)  (n-2) 

n 

I  y.[nmYZ-y.z.]z.-I[-y.]y.z.2 


1=1 


1=1 


n(n-l)  (n-2) 


n  n  n 

■nmV7  I  y.z.  +  I  y2z2  +  I  v2z2 
Yz  i=iJl  1     i=i Ji  1     j=i  'i  ) 


n(n-l)  (n-2) 


-nmYZ  nmyz  +  nmY2z2  +  nmY2z2 


n(n-l)  (n-2) 


2mY2z2  -  nmyz 


(n-1)  (n-2) 

E[Y.Z.YmZ.]  andE[Y.Z 
derived  similar  to  Eq.  17 


[17] 


j    m  mJ 


n    n  n 


E[Y.Z.Y  Z]  =hL  J=i  Y,W  -  2m^2 '  nm^z 
1  '  m  1         i^m   (n-1)  (n-2) 


[18] 


n(n-l)  (n-2) 


and 


n    n  n 

y  y  y  y  y  z.z 

E[Y.ZY  Z  ]  =  i=l  =1  m=l 

n(n-l)  (n-2) 


i    ]    m  mJ 


-  2mY2Z2  '  nmYZ 

(n-1)  (n-2) 


[19] 


Finally,  E[Y.Z.YmZk],ioi  m*k,  i*m,  i±k,j*m,  and 
j*k,  is  derived  as  follows: 


n    n    n  n 


I  I  I  I  yiymzizk 


E[YZ.Y  ZJ  =  i=i  Pim=ik=i 

1    1    m    k  i*j*m*k 


n(n-l)  (n-2)  (n-3) 


n    n  n 

J  1 1  ^ymZj  [-zfzf z  J 

i=l  j=l  m=l 
i*j*m 

n(n-l)  (n-2)  (n-3) 


nnn  nnn  nnn 

_-I 1 1  y^z,  -III  y^z2  -  II I  yiymzjZm 

—  i=l  )=1  m=l  i=l  j=l  m=l  i=l  j=l  m=l 

i*j*m  i*j*m  i*j*m 


n(n-l)  (n-2)  (n-3) 


[20] 


From  Eqs.  18  and  19,  1 1 1  =111  y,ymzizr 

1  i=l  i=l  m=l  i  =  l  j  =  l  m=l 


i*j*m 


i*j*m 


in  Eq.  20,  and 
E[Y.Z.Y  ZJ 

i    1    m  kJ 


nnn  nnn 

•2.III  7iW  -  III  yyj* 

i=l  j=l  m=l  i=l  j  =  l  m=l  1 
:  i*j*m  i*j*m  

n(n-l)  (n-2)  (n-3) 


[21] 


Using  the  results  of  Eqs.  16  and  18  and  substituting 
in  Eq.  21  we  get: 


E[YZY  ZJ  = 

i    i     m  kJ 


are 

2mY2Z2-nmYZ 

2mY2Z2-nmY2  mz2 

-2 

(n-1)  (n-2) 

n(n-l)(n-2) 

(n-1)  (n-2) 

n(n-l)(n-2) 

n(n-l)  (n-2)  (n-3) 

-2[2mY2Z2  -nmYZ]  -  [2mY2z2  -  nmy2  mz2] 
(n-1)  (n-2)  (n-3) 


2nmyz  +  nmY2  mz2  -  6mY2z2 


(n-1)  (n-2)  (n-3) 


[22] 


First  Two  Moments  of  Moran's  Bivariate  7YZ  Statistic 

The  first  two  moments  of  Jyz  may  be  determined 
under  the  null  hypothesis  of  no  spatial  autocorrelation, 
conditional  upon  the  n  observed  locations,  and  using 
the  expected  values  that  were  derived  in  the  previous 
section.  The  expected  value  of  7YZ  in  Eq.  1  is  determined 
in  Eqs.  23  to  25: 


E[/YZ]  = 


n  n 

I  I  w  .  Y.Z. 

i  =  l  j=l 


WVmY2  mz2 


[23] 
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Since  the  variables  w..,  W,  mY2,  and  mz2  (Eqs.  1,  7,  and  8) 
are  known  constants ,  or  invariant  under  random  permu- 
tations in  Eq.  23: 


n  n 


i=l  ,=1 

E[/YZ]  =  _^ 


I  I  w.EIY.Z,] 


\\f^mY2  mz2 


[24] 


The  E[  y.Z]  in  Eq.  24  has  the  same  value  for  all  pairs 
of  locations  i  and  ;',  and  using  the  definition  of  Win  Eq. 
1  and  the  value  of  E[y.Z]  derived  in  Eq.  12: 


E[/YZ]=li=i  |-i 


m 


YZ 


n-1 


\YVmY2  mz2 


w 


m. 


n-1 


YZ  = 


-  m 


YZ 


-  -P 


YZ 


\YVmY2  mz2 


(n-1)  VmY2  mz2 


n-1 


[25] 


where  mYZ,  my2,  and  mz2  are  defined  in  Eqs.  6,7,  and  8, 
respectively,  and  pYZ  is  a  measure  of  the  linear  correla- 
tion between  the  variables  Y  and  Z. 

E[/y ,  which  is  needed  to  compute  the  variance 
of  7YZ,  is  derived  from  Eq.  1  as  follows: 


E[/y  =  e 


=  E 


m=lk=l 

\YVmy2  rnz2 


n   n       n  n 


I    I      I    I       W«W«k  E[YiZiYmZJ 
i=l  j=l  m=lk=l 
(i*j)  (m*k) 

W2mY2  mz2 
The  numerator  of  Eq.  26  equals: 

i=l  j=l  m=lk=l 
(i*j)  (m*k) 


[26] 


n   n      n  n 

i*m,  j*k 
j*m,  i*k 


n  n  n  n 
I  I    I  I 

i*m.  j*k 
(j=m),  i*k 


n   n       n  n 

I  I     I  I 

+  i=l  i=l    m=lk=l  W..W  ,E[Y.Z.YZV] 
{*j,m*k  i)    mk         1    j    m  k 

i*m,  j*k 
j*m,  (i=k) 

n   n      n  n 

X  X     X  X 

+  i=l  i=l  m=lk=l  w..w  Jl[Y.Z.Y  ZJ 

i  .      \i  li    mk    L    i    i    m  kJ 

i*j,  m*k  ' 
(i=m),  j*k 
j*m,  i*k 

n  n      n  n 

X  X  X  X 

+ 1=1  i=i  m=ik=i  w..w  ,E[y.z.y  zj 

II  J  .A  m~J-K  1       li    mk    L    i    ]    m  kJ 
i*j,  m*k 

i*m,  (j=k) 
j*m,  i*k 

n   n      n  n 

X  X  X  X 

+  i=l  i=l  m=lk=l   W..W  ,E[Y.Z.Y  ZJ 
1   'J.,  i]    mk    L    l    )    m  kJ 

m*k 

i*m,  i*k 

(j=m),  (i=k) 

n  n      n  n 
I  I    I  I 

+  i=ij=i  m=ik=i  w  wmkE[YZYmZk] 

i*i,  m*k  1  1 

(i=m),  (j=k) 
j*m,  i*k 


[27] 


Since  any  other  summation  over  i,  j,  m,  and  k  in  Eq.  27 
is  infeasible  for  htj  and  m*k,  Eq.  27  simplifies  to  the 
following: 


x  x  x  x  ^mkE[rizjymzJ 

i=l  j=l  m=lk=l 
(i*j)  (m*k) 

=  I  i  £  i  ^^EtY^^ZJ 

i=l  j=l  m=l  k=l 
i*]>m*k 

+  I  I  I  ^"^[YZYZJ 

i=l  j=l  k=l 
i*j*k 

+  y  yy  ww  E[Y.Z.Y  Z] 

Z-    Z*    Z-        i]     mi    1     i    ]     m  iJ 

i=l  j=l  m=l 
l^m 

+IXX  ^^kEtWJ 

i=l  j=l  k=l 
i*j*k 

+  y  y  y  w..w  E[r.r  Z2] 

Zw    Zw    Zw         ]  i     mi  i     m  / 

i=l  j  =  l  m=l 


+  y  y  w..w.E[Y.Z.Y.Z.] 

Z-   Z-        ij    )i    1    i    ]    )  iJ 


+  y  y  ^..w.Etra2] . 

Z*    Z*        ij     ij    1     i    ;  j 

i=i if1 


[28] 


Next,  the  expectations  in  Eqs.  13-18,  and  22  can  be 
substituted  into  Eq.  28  to  produce  the  following: 
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X  X   X  X   w..w  .  E[Y.Z.Y  ZJ  = 

i=l  j=l  m=lk=l 
(m*k) 


n    n       n  n 


IIS  Iwiiwmk  2nm'z  +  nmTimzi-6mA; 
<i=1|S^k=1        /       (n-1)  (n-2)  (n-3) 


Computational  Simplifications 

The  formula  for  the  E[/y  in  Eq.  30,  which  is  needed 
to  compute  the  variance  of  Iyz,  can  be  further  simplified, 
thus  greatly  reducing  the  number  of  computations. 
First,  simplifying  notation  is  introduced,  and  then  used 
in  Eq.  30.  The  notation  in  Eqs.  31  to  34,  44  is  adopted 
from  Cliff  and  Ord  (1973,  pp.  25-26): 


n   n  n 


+  11  lKjwjk  +  wijwJ 


i=i j=i  k=i 

i*j*k 


+  (£  X  X  Cw..wik  +  W..W k.) 


i=l j=l  k=l 
i*j*k 


2mY2z2-nm^ 
(n-1)  (n-2) 


2mY2z2-nmY2mz2 
(n-1)  (n-2) 


n  n 

+  /  v  v  w;:wj     nmYZ  -  my2z2 


X  XVji] 
i=i  if1 


n-1 


+   X  X  w2.  nmY2m22-mY2z2 


[291 


i=i  j=i 


n-1 


Finally,  Eqs.  26  and  29  can  be  combined  to  compute 


E[iy 


n   n      n  n 


IX  X  X   X  wijwmk  1  |n[2mY2+  nyny]-  6mY2z2] 

li=l  j=l  m=l  k=l 
i*j*m*k 


(n-1)  (n-2)  (n-3) 


W2mY2mz2 


w..  -  yw.. 

i  l) 


W..  -  YW.. 
]      ^  i) 

i=l 


11         11  " 

S  =X  X   w.{w..  +  w.)  =  X  X  [w2.  +  W..W.) 

1  tltl   "  "     11    i=l  1=1   1     11  " 


[31] 
[32] 
[33] 


S2  =  XK-  +  w..)2  =  £(wf  +  2  w.w..  +  w.2.)  .  [34] 

i=l  i=l 

In  addition,  we  introduce  the  following  notation: 


S  = 


n  n 

y  y  w.  w.. 
i=l j=l 


s4=  X  Xw; 


i=l j=l 


y  w.w.. 

i=l  1  1 


[35] 


[36] 


[37] 


+  f£  X  X   KjWjk  +  WijWJ    ]  (2mY2z2-nmYZ 

(n-1)  (n-2) 


li=i  j=i  k=i 

i*j*k 


W2mY2mZ2 


n    n  n 


+  X  X  X   KiWik  +  WiiWJ  2mY2z2-nmY2mz2 


i=i  j=i  k=i 
i*j*k 


(n-1)  (n-2) 


W2mY2mz2 


n  n 

+  /Y  V  w.w  1  /nm2  -mv2,2 


X  X   "ii"H|  I  rz  — rz° 

i=1J=1         /  \      n-1       /  \W2mY2mz2 


+    X  X 


wz 


nmY2mz2  -  mY2z2| 


i=i  j=i 


n-1 


W2mY2mz2 


[30] 


S6=  X(w?  +  w2) 


[38] 


i=l 


It  can  be  verified  algebraically  from  Eqs.  35-38  that  S1 
and  S2  in  Eqs.  33  and  34,  44  equal: 


S,  =  2S  +  S„ 


[39] 
[40] 


We  have  found  it  useful  to  implement  both  the 
univariate  and  bivariate  Moran's  J  statistics  using  a 
computer  matrix  language.  Therefore,  the  computa- 
tional formulae  for  W,  Sv  S2,  S3,  S4,  S5  and  S6  are  given 
in  matrix  notation: 


w  = 

l'Wl, 

[41] 

s3  = 

l'(W<g>W')l, 

[42] 

s<  = 

l'(W®W)l, 

[43] 

s5  = 

(Wl)'(W'l)  = 

:  l'W  W'l  =  l'WWl, 

[44] 

s&  = 

(Wl)'(Wl)  + 

(W'l)'(W'l)  =  l'WWl 

+  l'WW'l  = 

l'(W'W  +  WW')l, 

[45] 
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51  =  S4  +  S3  =  l'(W®W+  W®  W')l,  [46] 

52  =  2S5  +  S6  =  l'(2WW  +  WW  +  WW')1,  [47] 
where 

1   =     is  a  nxl  vector  in  which  all  elements  are 
equal  to  1 

W  =    is  the  J7XI7  matrix  in  which  the  i/th  ele- 
ment equals  w..  and  all  diagonal  ele- 
ments equal  zero  (i.e.,  w..  =  0  for  all  i) 
W  =    is  the  matrix  transpose  of  W 
A®  B  =  denotes  element  by  element  multiplica- 
tion (i.e.,  the  7/th  element  of  A®  B  is 

equal  to  a.  b..). 

i]  ij' 

The  equality  in  Eq.  46  is  derived  from  Eqs.  33,  35, 
and  36,  and  the  equality  in  Eq.  47  is  derived  from  Eqs. 
34,  44,  37  and  38. 

The  first  term  in  Eq.  30  can  be  simplified  to  the 
following  using  the  definition  of  W in  Eq.  1,  and  Eqs.  31 
and  32: 


n   n     n  n 


n  n 


I  Z  I  I  W..W  . 
U=lj=lm=lk=l    11  mk 


/    n  n 

=IIw       I I  V 

i=l  j=i        I  m=l  k=l 
i*j  \  m*k*i*j 


mk 


n  n 
I  IW 


W 


[48] 

where  w..  is  an  element  of  the  ith  row  and  /th  column  of 
the  W  matrix  [W  =  l"Wl  Eqs.  41  to  45).  In  evaluating  Eq. 
48,  w..  and  w..  must  be  subtracted  only  once  from  W 
in  Eq.  48.  This  leads  to  the  following  result: 


n   n     n  n 

I  I  I  I  lw..w 

li=l  j=l  m=l  k=l 
i*j*m*k 


ij  mk 


n  n 
W2  -  I   (w..  +  W..)  W..  -   I  (w..  +  w..)w..  +  S, 
1=1  j=l 


n  n 
W2  -  I    (w.2.  +  W..W.)  -    I  (w..w..  +  W.2)  +  S„ 
i=l  i=l 


=  W2  -  I  (w2.  +  2w..w..  +  w..)  +  S, 


1=1 


=  w  -  s2  +  s, 


[50] 


This  agrees  with  the  derivation  by  Cliff  and  Ord  (1973, 
p.  26,  Eq.  2.18). 

The  next  term  in  Eq.  30  is  derived  as  follows: 

III  +  w^J I  I  w.iK  +  wj 


ti=l  j=l  k=l 

i*j*k 


ii=l  j=l  k=l 
i*j*k 


n  n 

n 

I   I  Wij 

I 

i=l j=l 

k=l 

\  k*i 

k*i 

M 

[si; 


Using  the  notation  in  Eqs.  31  and  32,  Eq.  51  simplifies 
to: 

i  i  i  Kwik +  wijwk1) 
,i=i  j=i  k=i 


n  n 

=  y  y  w..[(w.  -w..-w..)  +  (w..-w..-w..)] 

i=i if1 
1*) 


[52] 


Since  w..  =  w..  =  0  by  definition  in  Eq.  1,  Eq.  52  reduces 
to: 


n  n 


£  XW..[W  -  (w..)  -  (w..)  -  (w..-w..)  -(w..-w..)]  [49] 


where  w..  and  w..  in  Eq.  49  are  the  sums  of  the  ith  and  /th 
rows  of  the  W,  and  w..  and  w..  are  the  sums  of  the  ith  and 
/th  columns  of  W.  Next  if  we  use  the  definition  of  Win 
Eq.  1,  Sj  in  Eq.  33,  and  S2  in  Eq.  34  we  get: 


n   n  n 
III    (w  w   +  W  W,.) 
i=l  i=l  k=l  ' 

i*j*k 


n  n  n  n 


n         n         n  n 


Iw   Iw  -  I  Iww    +  Iw..Iw-I  Sww 

j=l    1    i=l         i=lj=l          '/     \i=l      ]=1       i=l  )=1  I 


J 


n        /  n   n     n  n 


I  I  wij  w  - 1 K- +  w-i)  iwij  ill  iwi 


w 


1=1  ]=1 
1*J 


i=l 


j=l      li=i  j=i  m=l  k=l 
i*j*m*k 


mk 


n  n  n  n 

Iw..(w..-w..)  +  I  w..(w..-w..)  -  2  11  w..w.. 

j=l  i=l  i=l j=l 


n  n  n  n 

=    -I  (w..+w.  )  Iw..  +  I  Iw..(w..+w  .) 


n  n 


2  y  w  .w..  -  2  I  I  w..w.. 
•^i     i  ij  ji 


[53] 


i=l 


i=l  j=l 
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From  Eqs.  35  and  37,  Eq.  53  simplifies  to  the  following 
term  in  Eq.  30: 


n   n  n 

III  (w  w  +  w  wj  I  =  2S5  -  2S3 

i=l  =1  k=l  J 
i*j*k 


[54] 


The  next  term  in  Eq.  30  is  derived  similar  to  Eqs.  51-  54 
using  Eqs.  36  and  38  as  follows: 


iii  Kjwik  +  wijwkj)  H  x  i  i  w  (wik+w ) 

i=l  j=i  k=l  J     \i=l  j=l  k=l 

i*j*k  /      \  i*j*k 


I  w.k  +    I  W 


11  11 

=  T,  I  w..[(w..-w..-w..)  +  (w.-w.-w..)] 


i=i  j=i 


n  n  n   n        \      /  n         n  nn 

=  [Iw..£w..-S  Iwf.  +  Iw..Iw.. -I  Iw 

1  i=l        j=l          i=l  j=l      /     \  i=l      i=l       i=l  j=l 


y w  .fw..-  w.  i  +y  w..(w..-  w.) -  2  y  y  w2. 

i=l  j=l  i=l j=l 


Using  S3  +  S4  =  S1  from  Eq.  46  and  2S5  +  S6  =  S2  from  Eq. 
47  and  substituting  in  Eq.  57,  and  rearranging  terms  we 
get: 

E[/y  = 


m2zn[2(W2-S2+S1)  -2(S5-S3)(n-3)  +S3(n-2)(n-3)]  \ 
-  m^MW-S^SJ  -2(S5-2SX)  (n-3)  +S1(n-2)  (n-3)] 
v+  my2mz2n[  (W2-S2+SJ  -  (S6-2S4)  (n-3)  +  S4(n-2)(n-3)]  / 


(n-1)  (n-2)  (n-3)  W2my2mz, 


[58] 


The  variance  of  IYZ  is  defined  as: 


Var(/YZ)  =  E[/y  -  E[/YZ]2  .  [59] 

E[Jyz]  in  Eq.  59  is  already  derived  in  Eq.  24,  and  E[/y  is 
expressed  in  Eq.  58,  therefore: 


/ 


+ 


Var(/YZ)  = 


mkil_V2(W2-S2+S1)  +  (2S3-2S5)(n-3)  +S3(n-2)(n-3j]  ^ 


my2mz2 


-m 


Y2Z2 


.  mY2mz2  / 


[6(W2-S2+S1)+  (4Sr2S2)  (n-3)+S1(n-2)  (n-3)] 


\        +  n  [  (W2-S2+S1)  +  (2S4-S6)  (n-3)  +  S4(n-2)(n-3)]  j 
(n-1)  (n-2)  (n-3)  W2 


n  n        \  n  n 

=  [  I  w2.  +1  w2   -  2  I  I  w? 

L=l  j=l        J  i=lj=l 


=     -  2S.  . 


[55] 


E[iy 


V 


Substituting  Eqs.  35,  36,  50,  54,  and  55  for  E[/y  in  Eq. 
30: 

/(W2-  S2+SJ  (2nm2z-  6my2z2+  nmy2mz2)\ 
+  (2S5-2S3)  (2mY2z2  -  nm2z)  (n-3) 
+  (S6  -2S4)  (2mY2z2  -  nmy2mz2)  (n-3) 
+  S3  (nm2z  -  my2z2)  (n-2)  (n-3) 
+  S4  (nmY2  mz2  -  my2z2)  (n-2)  (n-3)  J 
(n-1)  (n-2)  (n-3)  W2mY2mz2  [56] 
Collecting  terms  for  myz,  mY2z2,  and  mY2mz2  in  Eq.  56: 

Etiy  = 

mYZ[(W2-S2+S1)2n  -(2S5-2S3)n(n-3)  +S3n(n-2)(n-3)] 
+  mY2z2[-6(W2-S2+S1)  +2(2S5+S6-2S3-2S4)  (n-3)  -(S3+S4)  (n-2)  (n-3)]] 
h-  my2mz2[( W2-S2+S1)  n  -  (S6-2S4)  n(n-3)  +  S4n(n-2)(n-3)] 


(n-1)  (n-2)  (n-3)  W2my2mz, 


[57] 


m2 

YZ 


in-iy 


[60] 


In  most  applications,  the  spatial  weights  are  sym- 
metric, i.e.,  w..  =  w..  and  W  =  W.  In  this  special  case: 

Sa  =  2  l'(WxW)l 


52  =  4  l'(WW)l 

53  =  S4  =  SJ2 

55  =  S2/4 

56  =  S2/2 

Var(/yz)  = 


/  /mYZn 


2(W2-S2+S1)  +  (Sr^)(n-3)  +S1(n-2)(n-3) 


l  mY2mz2  , 

+  [5^rN)[6(W2-S2+S1)+  (4Sr2S2)  ^-3)45,(11-2)  (n-3)] 


\ 


+  n 


(W2-S2+Sa)  +  (St--J)  (n-3)  +  Sa(n-2)(n-3) 


(n-1)  (n-2)  (n-3)  W2 
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'  mYZ 

imy2mz2 


(n-1)2 


[61] 
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Univariate  Moran's  /  as  a  Special  Case 

Moran's  bivariate  JYZ  should  include  Moran's 
univariate  /  as  a  special  case,  as  presented  by  Cliff  and 
Ord  (1973,  pp  32-33).  The  purpose  of  this  section  is  to 
support  this  supposition  by  demonstrating  that  the 
mean  and  variance  of  Moran's  bivariate  JYZ  is  identical 
to  those  of  Moran's  univariate  /  when  the  two  response 
variables  Y  and  Zare  identical  (e.g.  Y=  Z). 

In  the  univariate  case,  y  =  z  for  all  locations, 

j  v  P 

1  <  p  <  n.  Substituting  this  identity  into  Eqs.  6,  7, 
8,  and  9: 


rn    n  n 


n    n  n 


+  I  I  I  w  w.k+X  I  Iw..wkil  (2m4-nm2m2 

\i=l  j=i  k=l  i=l  j=l  k=l  /  \7  7T7  7y\ 

L]*k  i*j*k  /  \(n-l)(n-2) 


[67] 


_  x 

mYZ  ~  P=l  ZpZp 

n 


_     I  ~2 

n 


=  Var(Z)  =  m 


[62] 


By  switching  the  i  and  /  subscripts  in  Eq.  67  for 

n    n     n  n    n  n 

III  w  w  and  I  I  X  w  w  we  get: 
i=ij=ik=i  i=i  j=i  k=i 

i*j*k  i*j*k 


my2 


n 

^  Z2 

n 


=  Var(Z)  = 


[63] 


mzz 


n 

I  Z2 
P=l  P 

n 


=  Var(Z)  =  m 


[64] 


n    n     n  n 


=  1111  w,wmk|  3nm2-6m4 


li=l  j=l  m=l  k=l 
i*j*m*k 


(n-1)  (n-2)  (n-3) 


mY2Z2 


^L,Z2Z2  =Lz4 

p-1  p  p      p-1  p 


=  m. 


n 


n 


[65] 


where  m2  and  m4  are  taken  from  the  notation  used  by 
CliffandOrd(1973,pp.  32-33).  Substituting  m2  and  m4 
(Eqs.  62-65)  in  Eq.  24  we  get: 


E['YZ]  = 


-m, 


-m„ 


(n-lWm,m,       (n-l)m,  (n-1) 


[66] 


2     2  v        '  2 

which  is  identical  to  the  results  of  Cliff  and  Ord  (1973, 
p.  32,  Eq.  2.31)  for  the  expected  value  of  Moran's 
univariate  J. 

Substituting  m2  and  m4  (Eqs.  62-65)  in  Eq.  29: 


E 


n   n      n  n 


=  11    II  w.  w  tE[Y.Z.Y  ZJ 

i=lj=lm=lk=l    11    mk       1   1    m  k 


+ 


in    n     n  n    n     n  \ 

/III  wwik+I  I  Iwwkix 

i=lj=lk=l  i=lj=lk=l 
i*j*k  i*j?tk 


n    n  n 


n    n  n 


III  w  wik+I  I  lwj(wki 

i=lj=lk=l  i=ij=ik=l 
\     i*j*k  i*j*k 


+  11  w^+w..) 
li=i  if1 


and  then  combining  terms: 


n    n     n  n 

=  11  I  I  I  w  wmk 

>i=l  j=l  m=l  k=l 
i*j*m*k 


nm22-m4 
n-1 


2m4-nm22 
in-l)(n-2)/ 


3nm2  -  6m, 

2  4 

(n-1)  (n-2)  (n-3) 


n   n      n  n 


=  11    II  Wi^mk  n(2m2+m2m2)-6m4 


,i=l  j=l  m=l  k=l 
i*j*m*k 


(n-1)  (n-2)  (n-3) 


n    n  n 


n    n  n 


+  |I  I  I  w..w.k+  III  w..wk.  2m4-nm; 
^i=1fikk=1  i=1tjkk=1  /\(n-D(n-2) 


n    n  n 


+  fl  I  I  (W..+W..)  (wik+wki) 

i=l  j=l  k=l 
i*j*k 


2m4-nm2 
(n-l)(n-2) 


n  n 


I  iw^+w,) 


,i=1.]=1 


nm2  -  m4 
n^l 


[68] 
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which  is  identical  to  the  results  given  by  Cliff  and  Ord 
(1973,  p.  33,  Eq.  2.33)  for  E[/y . 

The  next  step  is  to  verify  that  Ef/^]  in  Eq.  58  agrees 
with  the  equation  for  the  same  quantity  in  Cliff  and  Ord 
(1973,  p.  33,  Eq.  2.34).  Substituting  m2  and  m4  (Eqs.  62- 
65)  into  Eq.  58: 

E[/y  = 

n[2(W2-S2+S1)  +  (2S3-2S5)(n-3)  +S3(n-2)(n-3)]  \ 
+  l^^-)[6(W2-S2+S1)+  (4Sr2S2)  (n-SHS^)  (n-3)] 

Illy2      Z2  ' 


+  n  [  (Y^-ySJ  +  (2S4-S6)  (n-3)  +  S4(n-2)(n-3)] 
(n-1)  (n-2)  (n-3)  W2 


n 


3(W2-S2+Sa)  +  (2S3+2S4-2S5-S6)(n-3) 
+  (S3+S4)  (n-2)(n-3) 


\ 


+ 


[6(W2-S2+S1)+  (4Sr2S2)  (n-3)+S1(n-2)  (n-3)] 


y  \mY2mz2y 


(n-1)  (n-2)  (n-3)  W2 


[69] 


Since  S3  +  S4  =  from  Eq.  46  and  2S5  +  S6  =  S2  from  Eq. 
47,  Eq.  69  simplifies  to: 


E[/y  = 


-m. 


i  + 


_Vz2_ 
lmY2mz2 


n[3(W2-S2+SJ  +  (2SrS2)(n-3)  +  Sa(n-2)  (n-3)] 
X|[6(W2-S2+S1)+  (4Sr2S2)  (n-3)+S1(n-2)  (n-3)] 


(n-1)  (n-2)  (n-3)  W2 


/ 


-m 


+ 


Y2Z2 


lmY2mz2 


n[3W2+  (n2-3n+3)S1-nS2] 
|[6W2+  (n2-n)  Sr2nS2] 


[70] 


(n-1)  (n-2)  (n-3)  W2 


Eq.  70  is  identical  to  Eq.  2.34  in  Cliff  and  Ord  (1973). 
This  demonstrates  that  the  variance  of  Moran's  bivariate 
7YZ  in  Eq.  60  equals  the  variance  of  Moran's  univariate  / 
under  the  assumptions  in  Eqs.  62-65. 


were  originally  used  to  detect  changes  in  net  basal  area 
growth  of  natural,  undisturbed  shortleaf  pine  stands  in 
Georgia  between  1961-72  and  1972-82.  Because  of  high 
disturbance  rates  the  number  of  sample  plots  used  to 
estimate  growth  varied  in  each  period  (Bechtold  et  al. 
1991);  127  plots  for  the  1961-72  period  and  40  plots  for 
the  1972-82  period. 

Data  used  in  the  analysis  included  the  natural  loga- 
rithm of  gross  annual  pine  basal  area  growth  per  acre 
(G),  site  index  (S)  (base  age  50),  natural  logarithm  of 
stand  age  (A)  (midpoint  of  10-year  class),  natural  loga- 
rithm of  the  number  of  trees  per  acre  (N),  ratio  of  pine 
basal  area  to  basal  area  of  all  species  (P),  and  natural 
logarithm  of  annual  basal  area  mortality  per  acre  (M). 
For  a  more  detailed  description  of  the  data  see  Bechtold 
et  al.  (1991). 

Simulation 

Permutation  procedures  were  used  to  estimate  the  mean, 
variance,  skewness,  and  P-value  of  Moran's  bivariate  I 


YZ 


under  the  null  hypothesis  of  no  spatial  autocorrelation 
between  two  variables.  Since  the  number  of  possible 
permutations  tends  to  be  large  for  even  small  data  sets, 
an  approximation  to  the  permutation  procedure  was 
used  to  obtain  estimates  of  the  population  parameters 
under  the  null  hypothesis.  This  was  accomplished  by 
randomly  assigning  the  response  variables  to  the  n 
distinct  geographical  locations  and  calculating  IYZ.  This 
process  was  repeated  200,000  times  with  replacement 
for  each  of  the  15  pairs  of  variables  associated  with  the 
two  growth  cycles.  In  calculating  Jyz  the  inverse  dis- 
tance between  plots  was  used  as  a  weighting  factor. 

The  simulations  were  used  to  detect  any  errors  in 
the  derivation  of  IYZ  (Eq.  25)  and  its  variance  (Eq.  60).  A 
percent  difference  for  IYZ  was  computed  as  the  differ- 
ence between  the  average  of  the  200,000  estimates  of  IYZ 
from  the  simulations  and  the  expected  value  (Eq.  25) 
under  the  null  hypothesis,  divided  by  the  expected 
value.  This  was  done  for  each  of  the  1 5  pairs  of  variables 
in  each  growth  period.  In  addition,  the  variance  of  the 
mean  (i.e.  simulation  variance)  was  also  computed  as 
the  variance  among  the  200,000  estimates  of  IYZ.  In 
general,  the  variance  of  the  mean  is  the  best  estimate  of 
the  variability  and  can  be  used  to  evaluate  whether  the 
derivation  of  the  variance  of  IYZ  given  by  Eq.  60  is 
correct. 


Statistical  Properties  of  Moran's 
Bivariate  /YZ  DATA 

The  statistical  properties  of  Moran's  bivariate  /YZ 
were  evaluated  using  a  subset  of  USDA  Forest  Service 
inventory  data  from  Bechtold  et  al.  (1991).  The  data 


Finally,  the  P-value  for  the  realized  value  of  IYZ  (say 
IQ)  was  calculated  as  the  proportion  of  the  7YZ's  from  the 
200,000  simulations  that  were  more  extreme  than  IQ,  or 
equivalently,  the  probability  under  the  null  hypothesis 
given  by  PUYZ  >  J0).  A  P-value  was  also  calculated  using 
the  normal  and  Pearson  type  III  distributions. 
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Comparison  of  the  Expected  Value  and  Simulation 
Mean 

The  expected  value  of  IYZ  (Eq.  25)  agreed  with  the 
simulation  mean  of  Moran's  bivariate  IYZ.  In  general,  as 
the  number  of  observations  increased  from  40  (second 
growth  cycle)  to  127  (first  growth  cycle)  the  magnitude 
of  the  percent  difference  in  IYZ  decreased  from  2.7 
percent  to  -0.10  percent.  The  largest  percent  difference 
was  observed  when  the  linear  correlation  between  any 
two  variables  was  near  zero  (e.g.,  pYZ~0).  Since  the 
expected  value  of  Moran's  bivariate  test  statistic  IYZ 
approaches  zero  as  the  linear  correlation  between  any 
two  variables  approaches  zero,  the  percent  difference 
tends  toward  infinity.  Thus,  this  trend  is  likely  a 
numerical  artifact  and  not  an  indication  of  an  error  in 
the  derivation  of  the  expected  value  of  Moran's  bivariate 
test  statistic. 

For  small  sample  sizes  [n  =  40)  the  simulation  mean 
tends  to  underestimate  the  expected  value  with  increas- 
ing magnitude  as  the  linear  correlation  between  any  two 
variables  increased  in  a  positive  direction.  This  trend 
was  not  significant,  however,  at  the  0.05  level  of  signifi- 
cance. As  the  sample  size  increased  [n  =  127)  the 
percent  difference  was  observed  to  be  independent  of 
the  linear  correlation  between  any  two  variables.  As  the 
sample  size  increases,  we  expect  the  percent  difference 
to  approach  zero. 

Ratio  of  Variances 

The  ratios  of  the  simulation  variance  to  the  variance 
of  Moran's  bivariate  Jyz  under  the  null  hypothesis  were 
all  less  than  1.  A  ratio  less  than  1  indicates  an  underes- 
timation of  the  true  variance,  while  a  ratio  greater  than 
1  would  indicate  an  overestimation.  As  the  sample  size 
increased  from  40  to  127  the  average  ratio  increased 
from  0.961  to  0.985,  respectively.  We  conjecture  that  as 
the  sample  size  increases  the  simulation  variance  will 
approach  the  true  value  given  by  Eq.  60. 

Skewness 

The  skewness  for  Moran's  bivariate  IYZ  is  given  in 
Table  1 .  Except  for  two  cases  (G-N  and  S-P)  the  absolute 
value  of  the  skewness  was  less  than  0.01  for  both  data 
sets  analyzed  in  this  study,  indicating  that  Moran's 
bivariate  IYZ  tends  to  be  normally  distributed  for  sample 
sizes  n  >  40.  Mielke  (1986)  points  out  that  if  the  abso- 
lute value  of  the  skewness  is  less  than  0.01,  one  can 
reliably  assume  a  normal  distribution  in  estimating  P- 
values.  In  the  case  when  the  absolute  value  of  the 
skewness  is  greater  than  0.01,  the  Pearson  type  III 
distribution  can  be  used  to  account  for  skewness  in 
calculating  P-values. 


Table  1 .  Skewness  of  Moran's  bivariate  /YZ  associated  with  testing  the 
null  hypothesis  of  no  spatial  autocorrelation  between  two  variables. 


Data  set  1 

Data  set  2 

Variables1' 

n=  127 

n  =  40 

G-S 

0.0026 

0.0062 

G-A 

-0.0037 

0.0004 

G-N 

0  0055 

0.01592 

G-P 

0.0046 

0.0082 

G-M 

-0.0014 

-0.0017 

S-A 

-0.0001 

0.0007 

S-N 

0.0010 

-0.0071 

S-P 

0.0028 

0.01192 

S-M 

-0.0007 

0.0084 

A-N 

-0.0013 

0.0034 

A-P 

-0.0013 

0.0030 

A-M 

0.0013 

-0.0018 

N-P 

0.0034 

0.0069 

N-M 

0.0047 

0.0088 

P-M 

0.0002 

0.0058 

See  the  text  for  definitions  of  the  variables. 


z  If  the  absolute  value  of  the  skewness  is  greater  than  0.01,  one 
cannot  reliably  assume  a  normal  distribution  in  estimating 
P-values  (Mielke  1986). 

P-value 

Table  2  compares  the  P-values  for  Moran's  bivariate 
IYZ  under  the  null  hypothesis  of  no  spatial  autocorrelation 
between  two  variables.  P-values  were  calculated  using 
the  normal  distribution,  Pearson's  type  III  distribution 
which  takes  into  consideration  the  amount  of  skewness 
(Table  1),  and  as  the  proportion  of  the  200,000  simulated 
IYZ's  that  were  more  extreme  than  the  realized  value  of 
IYZ.  Since  the  absolute  value  of  the  skewness  for  Moran's 
bivariate  L7  was  less  than  0.01  for  all  but  two  simula- 

YZ 

tions,  P-values  based  on  the  normal  and  Pearson  type  III 
distribution  did  not  differ  appreciably  from  one  an- 
other. Comparing  these  P-values  with  those  obtained 
from  the  simulation  indicate  that  both  the  normal  and 
Pearson's  type  III  distribution  provide  a  reasonable 
approximation  of  the  tail  probabilities  under  the  null 
hypothesis  of  no  spatial  autocorrelation. 

Number  of  Simulations  When  Using  Small  Sample 
Sizes 

The  results  of  this  study  demonstrate  the  asymp- 
totic normality  of  Moran's  bivariate  IYZ  under  the  null 
hypothesis  for  two  data  sets  with  sample  sizes  n  >  40.  In 
this  case,  it  is  sufficient  to  take  z=(/YZ-E[/YZ])/V(/YZ)1/2  as 
a  standard  normal  variate.  The  null  hypothesis  of  no 
spatial  autocorrelation  between  two  variables  defined 
by  Y.  and  Z.  is  rejected  if  I  z  I  >  za/2.  If  IYZ  is  significantly 
larger  than  E[JYZ],  this  would  indicate  a  positive  spatial 
autocorrelation  between  the  two  variables ,  while  values 
significantly  less  than  E[Jyz]  would  indicate  a  negative 
spatial  autocorrelation  between  the  two  variables. 
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Table  2.  Comparison  of  P-values  for  a  1 -tailed  test  of  the  null  hypoth- 
esis of  no  spatial  autocorrelation  between  two  variables. 


Variables1 


NormaP 


Simulation3 


Pearson4' 


Data  set  1 ,  n  =  1 27 


G-S 

.499 

.487 

.499 

G-A 

.313 

.317 

.311 

G-N 

.265 

.270 

.264 

G-P 

.293 

.273 

.292 

G-M 

.350 

.348 

.348 

S-A 

4.7x1 05 

7.0x10  5 

0.7x10 

S-N 

.204 

.198 

.203 

S-P 

.099 

.097 

.097 

S-M 

.411 

.404 

.412 

A-N 

.265 

.260 

.263 

A-P 

.089 

.083 

.088 

A-M 

.092 

.088 

.090 

N-P 

.439 

.421 

.439 

N-M 

.115 

.102 

.114 

P-M 

.156 

.149 

.155 

Data  set  2,  n  =  40 


G-S 

7.8x10 4 

30.X10"4 

6.1x10 4 

G-A 

.023 

.024 

.020 

G-N 

4.7x1 0-4 

30.X104 

.2x1 0"4 

G-P 

8.5x1 0-5 

90.X105 

7.2x10 5 

G-M 

.292 

.318 

.325 

S-A 

.052 

.046 

.048 

S-N 

.086 

.075 

.083 

S-P 

.100 

.098 

.098 

S-M 

.382 

.383 

.379 

A-N 

.075 

.065 

.071 

A-P 

.066 

.059 

.063 

A-M 

.317 

.298 

.312 

N-P 

.008 

.011 

.007 

N-M 

.116 

.098 

.110 

P-M 

.299 

.294 

.294 

v  See  the  text  for  the  definitions  of  the  variables. 

a  Calculated  using  the  expected  value  and  variance  of  Moran's 
Bivariate  lYZ  under  the  null  hypothesis  of  no  spatial 
autocorrelation. 

31  Based  on  200,000  Monte  Carlo  simulations. 

4/  Based  on  the  Pearson  type  III  distribution  using  the  mean, 
variance  and  skewness  from  a  Monte  Carlo  simulation. 


Cliff  and  Ord  (1981)  point  out  that  for  small  sample 
sizes  [n  <  50)  the  assumption  of  normality  may  not  be 
valid.  The  departure  from  normality  may  be  due  to 
several  factors:  1)  the  spatial  distribution  of  the  observa- 
tions within  the  study  area;  2)  the  weights  used  [w..];  3) 
the  distribution  of  the  variables  Y.  and  Z;  and  4)  the 
sample  size,  n  (Cliff  and  Ord  1981).  Thus,  for  small 
sample  sizes  we  recommend  using  a  Monte  Carlo  simu- 
lation to  approximate  the  empirical  distribution  of 
Moran's  bivariate  JYZ  under  the  null  hypothesis. 

To  estimate  the  mean  and  skewness  We  recommend 
at  least  150,000  simulations,  and  at  least  200,000  simu- 
lations for  estimating  the  variance.  This  is  based  on 
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Figure  1.  Relative  error  associated  with  estimating  a  range  of  P- 
values  under  the  null  hypothesis  of  no  spatial  autocorrelation  as 
a  function  of  the  number  of  simulations  and  a  sample  size 
n  =  40.  Figure  1a  depicts  the  P-values  0.0008,  0.008,  and  0.023; 
while  Figure  1b  depicts  the  P-values  0.052,  0.10,  and  0.30. 

evaluations  of  smaller  numbers  of  simulations  during 
the  progression  of  this  study.  In  spite  of  the  large 
number  of  simulations  required  to  estimate  the  mean, 
variance,  and  skewness  under  the  null  hypothesis,  stable 
estimates  of  the  P-value  can  be  obtained  using  a  smaller 
number  of  simulations.  Figure  1  gives  examples  of  the 
relative  error  associated  with  estimating  a  range  of  P- 
values  under  the  null  hypothesis  of  no  spatial 
autocorrelation  as  a  function  of  the  number  of  simula- 
tions and  a  sample  size  of  n  =  40.  For  P-values  greater 
than  0.05  reasonable  estimates  can  be  obtained  with 
5,000  simulations  while  P-values  less  than  0.05  require 
at  least  10,000  simulation  (not  shown). 


Summary 

In  this  paper,  we  derive  the  expected  value  and 
variance  of  a  bivariate  version  of  Moran's  J,  a  measure  of 
spatial  autocorrelation.  Results  from  a  Monte  Carlo 
study  indicate  that  for  moderate  to  large  sample  sizes  (n 
>  40)  Moran's  bivariate  Jyz  tends  to  be  normally  distrib- 
uted. The  assumption  of  normality  begins  to  deteriorate 
for  sample  sizes  less  than  40. 

For  small  sample  sizes  [n  <  40),  a  better  approxima- 
tion to  the  distribution  of  the  test  statistic  is  needed. 
One  possible  approach  is  to  derive  the  third  and  possi- 
bly the  fourth  moments  of  the  statistic.  This  has  not 
been  pursued  to  date  because  of  the  complexity  of 
evaluating  the  higher  order  moments  under  the  assump- 
tion of  randomization.  Until  the  problem  is  solved,  the 
use  of  Monte  Carlo  simulations  for  inference  with  small 
sample  sizes  would  be  preferable  to  that  of  the  normal 
approximation. 
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Moran's  J  statistic  has  been  used  by  ecologists  and  geographers 
alike  to  test  for  the  presence  of  spatial  autocorrelation  in  a  single 
variable  over  a  two-dimensional  plane.  In  this  paper,  we  provide  the 
derivation  of  the  expected  value  and  variance  of  a  bivariate  version 
of  Moran's  J  for  use  with  multivariate  data  under  the  assumption  of 
spatial  independence.  We  also  demonstrate  that  Moran's  univariate 
/  statistic  is  a  special  case  of  Moran's  bivariate  JYZ.  Results  of  an 
extensive  Monte  Carlo  study  show  that  the  expected  value  and 
variance  are  reliable  for  several  data  sets  with  moderate  sample  sizes 
[n  =  40  and  127)  and  varying  degrees  of  correlation  among  different 
bivariate  surfaces.  For  small  sample  sizes  [n  <  40)  at  least  5,000 
simulations  should  be  used  to  estimate  the  P-value  to  test  the  null 
hypothesis  of  no  spatial  autocorrelation  between  two  variables. 
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The  Rocky  Mountain  Station  is  one  of  eight 
regional  experiment  stations,  plus  the  Forest 
Products  Laboratory  and  the  Washington  Office 
Staff,  that  make  up  the  Forest  Service  research 
organization. 

RESEARCH  FOCUS 

Research  programs  at  the  Rocky  Mountain 
Station  are  coordinated  with  area  universities  and 
with  other  institutions.  Many  studies  are 
conducted  on  a  cooperative  basis  to  accelerate 
solutions  to  problems  involving  range,  water, 
wildlife  and  fish  habitat,  human  and  community 
development,  timber,  recreation,  protection,  and 
multiresource  evaluation. 


RESEARCH  LOCATIONS 

Research  Work  Units  of  the  Rocky  Mountain 
Station  are  operated  in  cooperation  with 
universities  in  the  following  cities: 


Albuquerque,  New  Mexico 

Flagstaff,  Arizona 

Fort  Collins,  Colorado* 

Laramie,  Wyoming 

Lincoln,  Nebraska 

Rapid  City,  South  Dakota 


'Station  Headquarters:  240  W.  Prospect  Rd.,  Fort  Collins,  CO  80526 


